AIM : To Visualize the satellite products (SST, KPAR) as well as the in situ sites in Lyttelton harbour.

KPAR is estimated using satellite measurements (MODIS AQUA) of absorption and backscattering coefficients. Their spatial resolution is 500m-pixel wide.
SST data come from MODIS AQUA with a nominal resolution of 1km but has been downscalled to 500m-pixel wide.
The in situ sites refer to cross-shore transects where ecological monitoring of presence/absence and abundance of macroalgae has been conducted by the Cawthron institute and related to the dredging of the harbour.

Document Note : The maps, plots and app within this document are interactive so make sure you give them a play like zooming in and out in the maps but also on the plots. Clicking on the legend allows to only select and display the time series needed.

KPAR

Mean
Availability of pixels

Mean
Availability of pixels
Conclusion and further Work

Scenario 1

KPAR
SST
Conclusion and further Work

Scenario 2

KPAR
SST

Scenario 1 vs Scenario 2

KPAR
SST

KPAR

The KPAR dataset consists of monthly means of the attenuation coefficient estimated using the absorption and backscattering coefficient using MODIS-AQUA measurements of the emergent flux, radiance leaving the water. The dataset runs from July 2002 to March 2019.

Mean

## Read site sheet
site <- read.csv(file=paste0(getwd(),'/Coords_sites.csv'))
##

## KPAR stack and mean
kpar_mean <- raster('A200207_201903_KPAR_Mean_Lyttelton_QAA.tif')
kpar <- stack('A200207_201903_KPAR_MO_Lyttelton_QAA.tif')
##

## Bathy of the zone NIWA grid
## Download Bathy, plot it with sites and transform as contour
bathy <- raster('Bathy_Lyttelton.tif')
bathy[bathy>0] <- NA
bathycontour <- rasterToContour(bathy)
bathycontour <- spTransform(bathycontour,crs("+proj=longlat +ellps=WGS84 +datum=WGS84 +towgs84=0,0,0"))
##

## Plot KPAR mean and sites on leaflet map
pal <- colorNumeric(rev(rainbow(10)), values(kpar_mean),
                    na.color = "transparent")

m <- leaflet(site) %>% setView(lng = 173, lat = -43.55, zoom = 10) %>%
  addTiles()  %>%# Print the map
  addScaleBar(position = "bottomright",options = scaleBarOptions(imperial=F)) %>% 
  addMouseCoordinates() %>%
  addMarkers(site, lat = ~lat,lng = ~lon, popup = ~name) %>%
  addRasterImage(kpar_mean,layerId = "Kd", col=pal, opacity = 0.8,project=T) %>%
  addLegend(pal = pal, values = values(kpar_mean),title = "Kd (m-1)") %>%
  addImageQuery(x=kpar_mean,layerId = "Kd",type='click',project=T,position='bottomleft') %>%
  addPolylines(data=bathycontour,color = "black",opacity=0.2,popup = bathycontour$level)
m

Figure 1 -Mean KPAR product around Lyttelton harbour, NZ.

Key observations:

Some high values for pixels close to shore. This over-estimation could be due to the ‘land adjacency effect’ which is due to the land being an near infra-red (NIR) bright target which create a NIR halo affecting the atmospheric correction (Matt explanation).
To think about a method to take different pixes into account for the estimation of the KPAR value a the site. Maybe a compromise to find between the distance to the site and the availability of the pixel. Also possible to take more pixels into account.

Availability of pixels

kpar_NA_sum <- raster('A200207_201903_KPAR_PixAvailability_Lyttelton.tif')

kpar_NA_sum_proj <- projectRaster(kpar_NA_sum,crs=crs("+proj=longlat +ellps=WGS84 +datum=WGS84 +towgs84=0,0,0"),method='ngb')

pal <- colorNumeric(rev(heat.colors(100)), values(kpar_NA_sum_proj),
                    na.color = "transparent")

m <- leaflet(site) %>% setView(lng = 173, lat = -43.55, zoom = 10) %>%
  addTiles()  %>%# Print the map
  addScaleBar(position = "bottomright",options = scaleBarOptions(imperial=F)) %>% 
  addMouseCoordinates() %>%
  addMarkers(site, lat = ~lat,lng = ~lon, popup = ~name) %>%
  addRasterImage(kpar_NA_sum_proj,layerId = "Number of pixels available", col=pal, opacity = 0.8) %>%
  addLegend(pal = pal, values = values(kpar_NA_sum_proj),title = "Number of pixels available") %>%
  addImageQuery(x=kpar_NA_sum_proj,layerId =  "Number of pixels available",type='click',position='bottomleft') %>%
  addPolylines(data=bathycontour,color = "black",opacity=0.2,popup = bathycontour$level)
m

Figure 2 - Availability of pixels for the KPAR product around Lyttelton harbour, NZ.

Key observations:

The “holes” offshore in the pixels are likely to be due because of the reprojection for leaflet. Nothing to worry about.
Fewer pixels available close to shore/sites. As said previous section, to think about a method that takes into account the availability of pixels and distance to site.

SST

The SST dataset runs from July 2002 to March 2019 and consists of monthly means. I don’t know from which satellite it is from but could retrieve that at the occasion.

Mean

## SST stack and mean 
sst <- stack('A200207_201903_SST_MO_Lyttelton.tif')
sst_mean <- raster('A200207_201903_SST_Mean_Lyttelton.tif')
##

## Plot SST mean and sites on leaflet map
pal <- colorNumeric(rev(rainbow(100)), values(sst_mean),
                    na.color = "transparent")

m <- leaflet(site) %>% setView(lng = 173, lat = -43.55, zoom = 10) %>%
  addTiles()  %>%# Print the map
  addScaleBar(position = "bottomright",options = scaleBarOptions(imperial=F)) %>% 
  addMouseCoordinates() %>%
  addMarkers(site, lat = ~lat,lng = ~lon, popup = ~name) %>%
  addRasterImage(sst_mean,layerId = "SST (degC)", col=pal, opacity = 0.8,project=T) %>%
  addLegend(pal = pal, values = values(sst_mean),title = "SST (degC)") %>%
  addImageQuery(x=sst_mean,layerId =  "SST (degC)",type='click',project=T,position='bottomleft') %>%
  addPolylines(data=bathycontour,color = "black",opacity=0.2,popup = bathycontour$level)
m

Figure 3 - Mean SST product around Lyttelton harbour, NZ.

Key observations:

Odd values in coastal pixels. Either too high? (in Lyttelton harbour and way out of bays) or too low? (within the bays).
As previously, to think about a method to get the SST estimated value by satellite to the sites.

Availability of pixels

sst_NA_sum <- raster('A200207_201903_SST_PixAvailability_Lyttelton.tif')

pal <- colorNumeric(rev(heat.colors(100)), values(sst_NA_sum),
                    na.color = "transparent")

m <- leaflet(site) %>% setView(lng = 173, lat = -43.55, zoom = 10) %>%
  addTiles()  %>%# Print the map
  addScaleBar(position = "bottomright",options = scaleBarOptions(imperial=F)) %>% 
  addMouseCoordinates() %>%
  addMarkers(site, lat = ~lat,lng = ~lon, popup = ~name) %>%
  addRasterImage(sst_NA_sum,layerId = "Number of pixels available", col=pal, opacity = 0.8,project=T) %>%
  addLegend(pal = pal, values = values(sst_NA_sum),title = "Number of pixels available") %>%
  addImageQuery(x=sst_NA_sum,layerId =  "Number of pixels available",type='click',position='bottomleft',project=T) %>%
  addPolylines(data=bathycontour,color = "black",opacity=0.2,popup = bathycontour$level)
m

Figure 4 - Availability of pixels for the SST product around Lyttelton harbour, NZ.

Key observations:

Less pixels available where dodgy values are suspected.

Conclusion and further Work

Because the sat data is available from 07/2002 to 03/2019, to define a time frame to refine the dataset. Use the surveys dates, maybe to focus on ecological data during the baseline period (19-20/01/17 BL1 to 5-8/12/17 BL3) and do another analysis with the data after dredging phase (DP1 only date to fell into sat data dataset range).
I could then regenerate monthly means within the new time frame and start exploring different scenarios for assigning a value to sites.

Different scenarios:

1 : Extract values (KPAR and SST) at site pixel -> Will be not realistic because of missing values due to shore and effect of land.
2 : Extract values at pixel wiht highest availability from cluster of 16 pixels adjacent to site location.
2’ : Extract values at closest pixel offshore above threshold of certain pixel availability.
3 : Extract values from cluster of pixels at certain distance of site and above threshold of pixel availability.

Scenario 1

In order to have estimates of KPAR and SST for each sites, we investigate here the scenario of extracting the time series of the pixel where the site is. This approach is supposed to fail as the sites are very coastal so prone to pixels unavailability due to cloud coverage as well as the land adjacency effect which leads to failure in the atmospheric correction process.

KPAR

## Work on Lat/Lon of site coordinates, reptroject to utm
# Reprojection of sites lon/lat into tmerc (EBED crs)
LatLong_site <- data.frame(Y=site$lat,X=site$lon)
names(LatLong_site) <- c("Y","X")
coordinates(LatLong_site) <- ~ X + Y # longitude first
proj4string(LatLong_site) <- crs("+proj=longlat +ellps=WGS84 +datum=WGS84 +towgs84=0,0,0")
utm_sites <- spTransform(LatLong_site, crs( "+proj=tmerc +lat_0=0 +lon_0=173 +k=0.9996 +x_0=1600000 +y_0=10000000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs"))
##

## Extract values of pixel availability at sites location
kpar_NA_sum_velox <- velox('A200207_201903_KPAR_PixAvailability_Lyttelton.tif')
kpar_NA_sum_sites <- kpar_NA_sum_velox$extract_points(utm_sites)
sst_NA_sum_velox <- velox('A200207_201903_SST_PixAvailability_Lyttelton.tif')
sst_NA_sum_sites <- sst_NA_sum_velox$extract_points(utm_sites)
##

## KPAR
kpar_velox <- velox('A200207_201903_KPAR_MO_Lyttelton_QAA.tif')
kpar_sites <- kpar_velox$extract_points(utm_sites)
kpar_sites_df <- data.frame(t(kpar_sites))

colnames(kpar_sites_df) <- site$name
rownames(kpar_sites_df) <- NULL
dateseq <- seq.Date(as.Date("2002/7/1"), by = "month", length.out = dim(kpar_sites)[2])
kpar_sites_df$Date <- dateseq

kpar_sites_tb <- as_tibble(kpar_sites_df)
kpar_sites_tb <- kpar_sites_tb %>% pivot_longer(-Date,names_to="Site",values_to='Kpar') %>% group_by(Date)

q <- ggplot(kpar_sites_tb,aes(Date,Kpar,color=Site)) + geom_point() + geom_line() + labs(y="Kpar (/m)", x = "Date")
ggplotly(q)

Figure 5 - Time series of KPAR (/m) at the 21 sites around Lyttelton harbour, NZ. The values are extracted from the pixels where the sites are (Scenario 1).

Relevant observations:

No values for sites: PL14, PL03, PB10, PB11, PB03, PB02.
Very frequent gaps in time series.
Most of the values are over 0.2 /m and under 1 /m. Possible overestimation.

## KPAR mean, sd, median, min and max values over period
summary_kpar <- kpar_sites_tb %>% group_by(Site) %>% summarise(Mean=mean(Kpar,na.rm=T),Sd=sd(Kpar,na.rm=T),Median=median(Kpar,na.rm=T),Min=min(Kpar,na.rm=T),Max=max(Kpar,na.rm=T))

## `summarise()` ungrouping output (override with `.groups` argument)

summary_kpar$PixAv <- kpar_NA_sum_sites
kable(summary_kpar)

Site	Mean	Sd	Median	Min	Max	PixAv
BP01	0.4807346	0.1439982	0.4621158	0.2263836	0.9213209	86
BP02	0.5245741	0.1420317	0.5049014	0.2587270	0.9719090	83
BP03	0.5477195	0.1388396	0.5400223	0.3038899	0.9552073	48
BP05	0.3739728	0.1198098	0.3444114	0.2331560	0.7454379	47
BP06	0.4466732	0.1221743	0.4181812	0.2577299	0.9109438	106
BP08	0.4077619	0.1305797	0.3795969	0.1555093	1.0791025	90
BP10	0.4126234	0.1417000	0.3839348	0.2120700	0.8721717	67
BP13	0.4603269	0.1148118	0.4427282	0.2227825	0.7604557	90
BP14	0.6887431	0.2607180	0.6435227	0.2820656	1.4770980	70
LH01	0.5052427	0.1413214	0.5029341	0.2142351	0.8423380	73
LH02	0.4657920	0.1432052	0.4895723	0.2473924	0.6401853	9
LH07	0.6919909	0.2042824	0.6676622	0.3572536	1.4077549	60
LH10	0.6233862	0.1094229	0.6374094	0.4760661	0.7800066	7
PB02	NaN	NA	NA	Inf	-Inf	0
PB03	NaN	NA	NA	Inf	-Inf	0
PB10	NaN	NA	NA	Inf	-Inf	0
PB11	NaN	NA	NA	Inf	-Inf	0
PL02	0.8797934	0.2362799	0.9601306	0.6019558	1.1715964	5
PL03	NaN	NA	NA	Inf	-Inf	0
PL14	NaN	NA	NA	Inf	-Inf	0
PL16	0.8977092	0.2799484	0.9074726	0.5237424	1.2508827	8

SST

sst_velox <- velox('A200207_201903_SST_MO_Lyttelton.tif')
sst_sites <- sst_velox$extract_points(utm_sites)
sst_sites_df <- data.frame(t(sst_sites))

colnames(sst_sites_df) <- site$name
rownames(sst_sites_df) <- NULL
dateseq <- seq.Date(as.Date("2002/7/1"), by = "month", length.out = dim(sst_sites)[2])
sst_sites_df$Date <- dateseq

sst_sites_tb <- as_tibble(sst_sites_df)
sst_sites_tb <- sst_sites_tb %>% pivot_longer(-Date,names_to="Site",values_to='SST') %>% group_by(Date)

p <- ggplot(sst_sites_tb,aes(Date,SST,color=Site)) + geom_point() + geom_line() + labs(y="SST (degC)", x = "Date")
ggplotly(p)

Figure 6 - Time series of SST (degC) at the 21 sites around Lyttelton harbour, NZ. The values are extracted from the pixels where the sites are (Scenario 1).

Relevant observations:

No values for sites: PL03, PB11, PB10, PB03, PB02.
Very frequent gaps in time series. Some sites only have a few points (PL14), whereas others have only a few missing values (All BP for instance).
Most of the values are over 10 degC (Winter) and under 20 degC (Summer). Seems relatively good, however let’s have a look at the mean for every stations and see how the gaps in time series will affect the mean.

## SST mean, sd, median, min and max values over period
summary_sst <- sst_sites_tb %>% group_by(Site) %>% summarise(Mean=mean(SST,na.rm=T),Sd=sd(SST,na.rm=T),Median=median(SST,na.rm=T),Min=min(SST,na.rm=T),Max=max(SST,na.rm=T))

## `summarise()` ungrouping output (override with `.groups` argument)

summary_sst$PixAv <- sst_NA_sum_sites
kable(summary_sst)

Site	Mean	Sd	Median	Min	Max	PixAv
BP01	16.72645	3.505492	16.85250	8.695	24.035	100
BP02	16.93592	3.504716	17.20750	9.870	23.340	102
BP03	16.67018	3.449685	17.54500	8.170	23.110	51
BP05	14.50368	3.624410	14.19000	8.715	19.910	19
BP06	16.46508	3.629748	17.31250	8.520	23.965	118
BP08	15.74805	3.667408	16.24000	7.805	22.180	107
BP10	16.96604	3.810822	17.72750	8.300	23.660	78
BP13	14.59495	4.071016	14.70000	7.100	22.625	77
BP14	16.18449	4.400136	16.31474	9.260	25.740	50
LH01	14.94987	4.823432	15.33750	7.460	25.455	74
LH02	13.81568	5.620958	12.02000	5.640	24.265	26
LH07	18.14572	4.455413	18.60000	8.575	26.540	79
LH10	22.14909	3.711519	22.14909	12.875	30.245	23
PB02	NaN	NA	NA	Inf	-Inf	0
PB03	NaN	NA	NA	Inf	-Inf	0
PB10	NaN	NA	NA	Inf	-Inf	0
PB11	NaN	NA	NA	Inf	-Inf	0
PL02	17.39167	4.701594	16.73083	9.650	24.750	10
PL03	NaN	NA	NA	Inf	-Inf	0
PL14	7.54500	0.000000	7.54500	7.545	7.545	2
PL16	18.87106	4.213544	20.63000	9.630	24.130	27

Conclusion and further Work

The relative high number of missing values as well as the possible impact of the adjacent land lead to the conclusion that the scenario 1-way of getting estimated values for each site does not seem to show realistic values.
Let’s carry on the exploration and analyse the other scenarios.

Scenario 2

In order to have estimates of KPAR and SST for each sites, we investigate here the scenario of extracting the time series from a pixel with the maximum availability (using the pixel availability raster) within a cluster of 16 pixels adjacent to the site location (17 pixels including site, see what’s the structure of the cluster).

KPAR

## KPAR
ngb_pix_maxPixAv_index_save <- c()
sites_kpar_mean_save <- c()
for (i in 1:dim(site)[1]) {
  sites_tmerc <- data.frame(utm_sites)[i,]
  sites_kpar_mean <- raster::extract(kpar_mean,sites_tmerc,cellnumbers=T) #Get cell number of pixel where site
  ngb_pixels <- raster::adjacent(kpar_mean,cells=sites_kpar_mean[,1],pairs=F,directions=16,include=T)
  
  ngb_pix_maxPixAv_index <- ngb_pixels[which(kpar_NA_sum[ngb_pixels]==max(kpar_NA_sum[ngb_pixels]))] #Return index of pixel in cluster of 16 pixels around that have max pixel availability
  # Issues with PB10 and PB3 -> Seems like no pixels available in vicinity
  if (length(ngb_pix_maxPixAv_index)>1) { 
    #Manually solve the issue affecting original index to pixel
    ngb_pix_maxPixAv_index <- sites_kpar_mean[,1]
  }
  ngb_pix_maxPixAv_index_save <- c(ngb_pix_maxPixAv_index_save,ngb_pix_maxPixAv_index)
}

# Position of the "new" pixels where extractions of scenario 2 is from
pix_index_maxPixAv <- xyFromCell(kpar_mean,cell=ngb_pix_maxPixAv_index_save,spatial=F)#Convert to sp feature, to use velox::
Tmerc_site_s2 <- data.frame(Y=pix_index_maxPixAv[,2],X=pix_index_maxPixAv[,1])
coordinates(Tmerc_site_s2) <- ~ X + Y # longitude first
proj4string(Tmerc_site_s2) <- crs("+proj=tmerc +lat_0=0 +lon_0=173 +k=0.9996 +x_0=1600000 +y_0=10000000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs")
LatLong_sites_s2 <- spTransform(Tmerc_site_s2,crs("+proj=longlat +ellps=WGS84 +datum=WGS84 +towgs84=0,0,0"))
LatLong_sites_s2$name <- paste0(site$name,"_s2")
LatLong_sites_s2_df <- data.frame(LatLong_sites_s2)

# Extract time series from "new" pixels
pix_index_maxPixAv_sp <- xyFromCell(kpar_mean,cell=ngb_pix_maxPixAv_index_save,spatial=T)#Convert to sp feature, to use velox::
kpar_sites_maxPixAv <- kpar_velox$extract_points(pix_index_maxPixAv_sp)#Classic extraction process
kpar_sites_maxPixAv_df <- data.frame(t(kpar_sites_maxPixAv))

colnames(kpar_sites_maxPixAv_df) <- LatLong_sites_s2$name#paste0(site$name[1],'_maxPixAv')
rownames(kpar_sites_maxPixAv_df) <- NULL
dateseq <- seq.Date(as.Date("2002/7/1"), by = "month", length.out = dim(kpar_sites_maxPixAv_df)[1])
kpar_sites_maxPixAv_df$Date <- dateseq

kpar_sites_maxPixAv_tb <- as_tibble(kpar_sites_maxPixAv_df)
kpar_sites_maxPixAv_tb <- kpar_sites_maxPixAv_tb %>% pivot_longer(-Date,names_to="Site",values_to='Kpar') %>% group_by(Date)

q <- ggplot(kpar_sites_maxPixAv_tb,aes(Date,Kpar,color=Site)) + geom_point() + geom_line() + labs(y="Kpar (/m)", x = "Date")
ggplotly(q)

Figure 7 - Time series of KPAR (/m) at the 21 sites around Lyttelton harbour, NZ. The values are extracted from the pixels with the highest availability in a cluster of 16 pixels adjacent to each sites. Scenario 2.

m <- leaflet(data=LatLong_sites_s2_df) %>% setView(lng = 173, lat = -43.55, zoom = 10) %>%
  addTiles()  %>%# Print the map
  addMarkers(LatLong_sites_s2_df, lat = ~Y,lng = ~X, popup = ~name)
m

Figure 8 - Position of the Pixels where time series are extracted in the scenario 2 case.

Relevant observations:

TO FILL (Less gap, Position of the new pixels, PB03 and PB10 no value coz no ngb)

## KPAR mean, sd, median, min and max values over period
summary_kpar_s2 <- kpar_sites_maxPixAv_tb %>% group_by(Site) %>% summarise(Mean=mean(Kpar,na.rm=T),Sd=sd(Kpar,na.rm=T),Median=median(Kpar,na.rm=T),Min=min(Kpar,na.rm=T),Max=max(Kpar,na.rm=T))

## `summarise()` ungrouping output (override with `.groups` argument)

summary_kpar_s2$PixAv <- kpar_NA_sum[ngb_pix_maxPixAv_index_save]
kable(summary_kpar_s2)

Site	Mean	Sd	Median	Min	Max	PixAv
BP01_s2	0.3736091	0.0757962	0.3536769	0.2344388	0.7235438	168
BP02_s2	0.3566422	0.0690696	0.3449551	0.2511407	0.7337249	184
BP03_s2	0.3382366	0.0698441	0.3276764	0.2026570	0.6407034	169
BP05_s2	0.3479043	0.0756010	0.3367585	0.2061284	0.6099731	163
BP06_s2	0.3381993	0.0727750	0.3218081	0.2134089	0.5809690	190
BP08_s2	0.3111466	0.0691485	0.2943282	0.1995939	0.6331935	190
BP10_s2	0.3616482	0.0944326	0.3516572	0.1610703	0.9613709	185
BP13_s2	0.3749535	0.0785294	0.3640774	0.1974660	0.6571044	180
BP14_s2	0.3566313	0.0928171	0.3355087	0.1544371	0.6814184	177
LH01_s2	0.4792453	0.1089673	0.4638782	0.2742959	0.8294032	138
LH02_s2	0.5040968	0.1202039	0.4734725	0.2543523	0.9888393	145
LH07_s2	0.5214943	0.1296387	0.4985805	0.2543523	0.9773747	146
LH10_s2	0.4840434	0.1156862	0.4757362	0.2742959	0.8221809	136
PB02_s2	0.4004347	0.0580121	0.4073745	0.3243970	0.4625929	4
PB03_s2	NaN	NA	NA	Inf	-Inf	0
PB10_s2	NaN	NA	NA	Inf	-Inf	0
PB11_s2	1.1407515	NA	1.1407515	1.1407515	1.1407515	1
PL02_s2	0.6564987	0.2188256	0.6817867	0.2699000	1.0174295	41
PL03_s2	1.0310847	0.2605226	1.1485419	0.6019558	1.3445294	11
PL14_s2	1.0310847	0.2605226	1.1485419	0.6019558	1.3445294	11
PL16_s2	0.4569534	0.1356378	0.4164062	0.2840720	0.9177421	87

SST

## SST
ngb_pix_maxPixAv_index_save <- c()
sites_sst_mean_save <- c()
for (i in 1:dim(site)[1]) {
  sites_tmerc <- data.frame(utm_sites)[i,]
  sites_sst_mean <- raster::extract(kpar_mean,sites_tmerc,cellnumbers=T) #Get cell number of pixel where site
  ngb_pixels <- raster::adjacent(sst_mean,cells=sites_sst_mean[,1],pairs=F,directions=16,include=T)
  
  ngb_pix_maxPixAv_index <- ngb_pixels[which(sst_NA_sum[ngb_pixels]==max(sst_NA_sum[ngb_pixels]))] #Return index of pixel in cluster of 16 pixels around that have max pixel availability
  # Issues with PB10 and PB3 -> Seems like no pixels available in vicinity
  if (length(ngb_pix_maxPixAv_index)>1) { 
    #Manually solve the issue affecting original index to pixel
    ngb_pix_maxPixAv_index <- sites_sst_mean[,1]
  }
  ngb_pix_maxPixAv_index_save <- c(ngb_pix_maxPixAv_index_save,ngb_pix_maxPixAv_index)
}

# Position of the "new" pixels where extractions of scenario 2 is from
pix_index_maxPixAv_sst <- xyFromCell(sst_mean,cell=ngb_pix_maxPixAv_index_save,spatial=F)#Convert to sp feature, to use velox::
Tmerc_site_s2_sst <- data.frame(Y=pix_index_maxPixAv_sst[,2],X=pix_index_maxPixAv_sst[,1])
coordinates(Tmerc_site_s2_sst) <- ~ X + Y # longitude first
proj4string(Tmerc_site_s2_sst) <- crs("+proj=tmerc +lat_0=0 +lon_0=173 +k=0.9996 +x_0=1600000 +y_0=10000000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs")
LatLong_sites_s2_sst <- spTransform(Tmerc_site_s2_sst,crs("+proj=longlat +ellps=WGS84 +datum=WGS84 +towgs84=0,0,0"))
LatLong_sites_s2_sst$name <- paste0(site$name,"_s2")
LatLong_sites_s2_sst_df <- data.frame(LatLong_sites_s2_sst)

# Extract time series from "new" pixels
pix_index_maxPixAv_sp <- xyFromCell(sst_mean,cell=ngb_pix_maxPixAv_index_save,spatial=T)#Convert to sp feature, to use velox::
sst_sites_maxPixAv <- sst_velox$extract_points(pix_index_maxPixAv_sp)#Classic extraction process
sst_sites_maxPixAv_df <- data.frame(t(sst_sites_maxPixAv))

colnames(sst_sites_maxPixAv_df) <- LatLong_sites_s2_sst$name#paste0(site$name[1],'_maxPixAv')
rownames(sst_sites_maxPixAv_df) <- NULL
dateseq <- seq.Date(as.Date("2002/7/1"), by = "month", length.out = dim(sst_sites_maxPixAv_df)[1])
sst_sites_maxPixAv_df$Date <- dateseq

sst_sites_maxPixAv_tb <- as_tibble(sst_sites_maxPixAv_df)
sst_sites_maxPixAv_tb <- sst_sites_maxPixAv_tb %>% pivot_longer(-Date,names_to="Site",values_to='SST') %>% group_by(Date)

q <- ggplot(sst_sites_maxPixAv_tb,aes(Date,SST,color=Site)) + geom_point() + geom_line() + labs(y="SST (degC)", x = "Date")
ggplotly(q)

Figure 9 - Time series of SST (degC) at the 21 sites around Lyttelton harbour, NZ. The values are extracted from the pixels with the highest availability in a cluster of 16 pixels adjacent to each sites. Scenario 2.

m <- leaflet(data=LatLong_sites_s2_sst_df) %>% setView(lng = 173, lat = -43.55, zoom = 10) %>%
  addTiles()  %>%
  addMarkers(LatLong_sites_s2_sst_df, lat = ~Y,lng = ~X, popup = ~name)
m

Figure 10 - Position of the pixels where time series of SST are extracted from in the scenario 2 case.

Relevant observations:

TO FILL (Less gap, Position of the new pixels, PB03 and PB10 no value coz no ngb, PL02 new pixel is in LH.)

## KPAR mean, sd, median, min and max values over period
# SST mean, sd, median, min and max values over period
summary_sst_s2 <- sst_sites_maxPixAv_tb %>% group_by(Site) %>% summarise(Mean=mean(SST,na.rm=T),Sd=sd(SST,na.rm=T),Median=median(SST,na.rm=T),Min=min(SST,na.rm=T),Max=max(SST,na.rm=T))

## `summarise()` ungrouping output (override with `.groups` argument)

summary_sst_s2$PixAv <- sst_NA_sum[ngb_pix_maxPixAv_index_save]
kable(summary_sst_s2)

Site	Mean	Sd	Median	Min	Max	PixAv
BP01_s2	14.55268	3.2726353	15.1200	7.785000	20.520	172
BP02_s2	16.93592	3.5047162	17.2075	9.870000	23.340	102
BP03_s2	14.35078	3.2165821	14.6975	8.120000	20.775	166
BP05_s2	14.69212	3.3178140	14.9900	8.479999	21.245	142
BP06_s2	13.86971	3.3077064	13.9600	8.264999	19.750	194
BP08_s2	13.94149	3.2028660	14.2700	8.160000	19.295	187
BP10_s2	13.93659	3.1511820	14.1350	8.255000	19.590	184
BP13_s2	14.47864	3.5653955	14.6300	8.264999	21.900	160
BP14_s2	14.44132	3.4180549	14.6200	8.040000	21.285	157
LH01_s2	14.78414	3.8983909	15.4800	7.535000	22.210	153
LH02_s2	14.93382	4.3012935	15.1625	7.610000	25.065	164
LH07_s2	14.82535	3.7880783	15.3200	7.580000	22.195	155
LH10_s2	14.95026	3.8885392	15.6750	7.770000	22.660	151
PB02_s2	13.69000	0.3400002	13.6900	13.349999	14.030	3
PB03_s2	NaN	NA	NA	Inf	-Inf	0
PB10_s2	NaN	NA	NA	Inf	-Inf	0
PB11_s2	21.83500	6.1904860	22.7275	13.520000	28.365	4
PL02_s2	19.21261	4.1895382	19.9050	9.525000	28.535	58
PL03_s2	19.00000	3.9163583	18.7625	13.145000	25.945	14
PL14_s2	17.37219	5.2478802	17.8000	5.345000	25.330	17
PL16_s2	16.41398	3.6061747	16.9075	8.139999	22.940	106

Scenario 1 vs Scenario 2

The following tables show the difference of the values from the scenario 1 minus scenario 2 for each sites. Negative values mean that the scenario 2 leads to an increase in the parameter/value from scenario 1.

KPAR

## KPAR - Relative difference: s1 - s2
is.na(summary_kpar) <- sapply(summary_kpar, is.infinite)#Replace INF value by NA
is.na(summary_kpar_s2) <- sapply(summary_kpar_s2, is.infinite)#Replace INF value by NA
rdiff_s1s2_kpar <- summary_kpar[,2:7] - summary_kpar_s2[,2:7]
rdiff_s1s2_kpar$Site <- summary_kpar$Site
kable(rdiff_s1s2_kpar)

Mean	Sd	Median	Min	Max	PixAv	Site
0.1071255	0.0682020	0.1084389	-0.0080553	0.1977770	-82	BP01
0.1679319	0.0729621	0.1599463	0.0075863	0.2381842	-101	BP02
0.2094828	0.0689955	0.2123459	0.1012329	0.3145038	-121	BP03
0.0260685	0.0442089	0.0076529	0.0270276	0.1354648	-116	BP05
0.1084739	0.0493994	0.0963730	0.0443210	0.3299748	-84	BP06
0.0966152	0.0614312	0.0852687	-0.0440846	0.4459090	-100	BP08
0.0509752	0.0472675	0.0322776	0.0509997	-0.0891992	-118	BP10
0.0853734	0.0362824	0.0786507	0.0253166	0.1033512	-90	BP13
0.3321118	0.1679009	0.3080139	0.1276285	0.7956796	-107	BP14
0.0259974	0.0323541	0.0390559	-0.0600609	0.0129349	-65	LH01
-0.0383048	0.0230013	0.0160998	-0.0069599	-0.3486540	-136	LH02
0.1704966	0.0746438	0.1690816	0.1029013	0.4303802	-86	LH07
0.1393428	-0.0062633	0.1616731	0.2017702	-0.0421743	-129	LH10
NaN	NA	NA	NA	NA	-4	PB02
NaN	NA	NA	NA	NA	0	PB03
NaN	NA	NA	NA	NA	0	PB10
NaN	NA	NA	NA	NA	-1	PB11
0.2232947	0.0174544	0.2783440	0.3320558	0.1541669	-36	PL02
NaN	NA	NA	NA	NA	-11	PL03
NaN	NA	NA	NA	NA	-11	PL14
0.4407557	0.1443106	0.4910663	0.2396705	0.3331407	-79	PL16

##

Key observations:

Considerable gain of pixel availability, except for 6 sites.
Not relevant methods for 6 sites where originally NA.
Scenario 2 leads to decrease of KPAR mean.

SST

## SST - Relative difference: s1 - s2
is.na(summary_sst) <- sapply(summary_sst, is.infinite)#Replace INF value by NA
is.na(summary_sst_s2) <- sapply(summary_sst_s2, is.infinite)#Replace INF value by NA
rdiff_s1s2_sst <- summary_sst[,2:7] - summary_sst_s2[,2:7]
rdiff_s1s2_sst$Site <- summary_sst$Site
kable(rdiff_s1s2_sst)

Mean	Sd	Median	Min	Max	PixAv	Site
2.1737668	0.2328568	1.7325006	0.9099998	3.5149994	-72	BP01
0.0000000	0.0000000	0.0000000	0.0000000	0.0000000	0	BP02
2.3194007	0.2331030	2.8475003	0.0500002	2.3349991	-115	BP03
-0.1884315	0.3065963	-0.8000002	0.2350006	-1.3349991	-123	BP05
2.5953728	0.3220413	3.3524995	0.2550001	4.2150002	-76	BP06
1.8065653	0.4645421	1.9700003	-0.3550000	2.8850002	-80	BP08
3.0294443	0.6596399	3.5924993	0.0450001	4.0699997	-106	BP10
0.1163033	0.5056204	0.0700002	-1.1649995	0.7250004	-83	BP13
1.7431729	0.9820810	1.6947432	1.2200003	4.4549999	-107	BP14
0.1657225	0.9250412	-0.1424999	-0.0749998	3.2450008	-79	LH01
-1.1181394	1.3196644	-3.1424999	-1.9699998	-0.7999992	-138	LH02
3.3203697	0.6673344	3.2800007	0.9949999	4.3449993	-76	LH07
7.1988307	-0.1770199	6.4740906	5.1050000	7.5849991	-128	LH10
NaN	NA	NA	NA	NA	-3	PB02
NaN	NA	NA	NA	NA	0	PB03
NaN	NA	NA	NA	NA	0	PB10
NaN	NA	NA	NA	NA	-4	PB11
-1.8209427	0.5120556	-3.1741657	0.1250000	-3.7849998	-48	PL02
NaN	NA	NA	NA	NA	-14	PL03
-9.8271873	-5.2478802	-10.2549996	2.1999998	-17.7850003	-15	PL14
2.4570789	0.6073695	3.7224998	1.4900007	1.1900005	-79	PL16

Key observations:

Considerable gain of pixel availability, except for 5 sites.
Not relevant methods for 5 sites where originally NA
Scenario 2 leads to decrease of SST mean except at 4 sites (BP05, LH02, PL02, PL14).
PB02 wtf, to check

Sat_Product_Viz

Francois Thoral

27 April 2020

Table of contents

KPAR

Mean

Availability of pixels

SST

Mean

Availability of pixels

Conclusion and further Work

Scenario 1

KPAR

SST

Conclusion and further Work

Scenario 2

KPAR

SST

Scenario 1 vs Scenario 2

KPAR

SST